% Esta funcion calcula los polos para una plantilla (La salida es el filtro NORMALIZADO! wp=1)
% Si le entra un n=0 significa que el usuario no especifico el n y la funcion calcula el optimo.
% Si le entra un n, se hace el filtro NORMALIZADO con ese n (Puede o no respetar la plantilla y el unico dato necesario sera el Ap)
% Ap y Aa estan en dB.
function [Z, P, gain, order,Q, wn] = filter_chevyshev_inv (filterType, aa, ap, wa, wp, waMENOS, waMAS, wpMENOS, wpMAS, order, Q)

%% calculo epsilon, n (si no fue espesificado) y a
epsilon=1/sqrt(10^(aa/10)-1);
     switch filterType
        case 2,
        wn=wp/wa;
        case 3,
        wn=wa/wp;
        case 1,
        wn=(wpMAS-wpMENOS)/(waMAS-waMENOS);
        case 4,
        wn=(waMAS-waMENOS)/(wpMAS-wpMENOS);
    end 

if order==0

    epsilon=1/sqrt(10^(aa/10)-1)
    order=ceil(acosh(1/(epsilon*sqrt(10^(ap/10)-1)))/acosh(1/wn))
end
if (Q~=0)
    Q2=0;
    order=0;
    while (Q2 <= Q)
    order=order+1;
    Q2=(1/2)*sqrt(1+(1/(tan(pi*(2*order-1)/(2*order))*tan(asinh(1/epsilon)/order)))^2);
    end 
    order=order-1;
end

a=asinh(1/epsilon)/order;

%% inicializacion variables
sigma=zeros(1,order);
w=zeros(1,order);
if (mod(order,2))
    Z=zeros(1,order-1);
else Z=zeros(1,order);
end
k2=1;


%% calcula los polos y ceros
for k=1:order
    titas=(2*k-1)*pi/(2*order);
    sigma(k)=-sin(titas)*sinh(a);
    w(k)=cos(titas)*cosh(a);
    
    if (k~=ceil(order/2) && mod(order,2))
        Z(k2)=1/cos(titas);
        k2=k2+1;
    elseif (~mod(order,2))
        Z(k2)=1/cos(titas);
        k2=k2+1;
    end
end

w0=sqrt(sigma.*sigma + w.*w);

alpha=sigma./(w0.*w0);
beta=-w./(w0.*w0);
Z=roundDecimals(Z);
Z=1i.*Z + 0.00001;   % Corrector para que la fase salga bien
P=alpha+1i*beta;

Q=abs(w0./(2*sigma));

% Calculo la constante
    gain=1;
    for i=1:length(Z)
        gain=gain/Z(i); 
    end
    for i=1:length(P)
        gain=gain*P(i);
        
    end
    gain=real(gain);


end